%~ \chapter{Aproximaci\'on de funciones}
\chapter{Cuadrados M\'inimos}

Se tienen ciertas mediciones con un determinado error y se conjetura 
que la fuente corresponde a una funci\'on determinada. 
Por ejemplo, se puede suponer que los puntos corresponden a una recta 
y se desea entonces hallar los valores de $c$ y $d$ para los cuales 
la recta $y = cx + d$ aproxima \textsl{mejor} a los puntos. 

La funci\'on $f$ se utiliza para aproximar a los puntos 
$p(x_0), \dots, p(x_m)$. 
Existen diversas formas de definir el concepto de 
\textsl{mejor aproximaci\'on}. 
\begin{itemize}
	\item \textbf{Desviaci\'on absoluta} $\displaystyle 
		\arg{\min{ \{ \sum_{i=1}^m |f(x_i) - p(x_i)| \} } }$ \\
		La funci\'on m\'odulo no es derivable en el cero y dificulta 
		los c\'alculos. 
		Adem\'as no necesariamente puede obtenerse una soluci\'on.
	\item \textbf{Mini-Max}	$\displaystyle 
		\arg{\min{ \{ \max_{1\le i \le m}|f(x_i) - p(x_i)| \} } }$ \\
		No puede obtenerse una soluci\'on mediante m\'etodos 
			elementales. 
		Otorga demasiada importancia a pocos elementos an\'omalos 
			(outliers).
	\item \textbf{Cuadrados M\'inimos} $\displaystyle 
		\arg{\min{ \{ \sum_{i=1}^m [f(x_i) - p(x_i)]^2 \} } }$ \\
		Esta forma facilita los c\'alculos.
\end{itemize}

\section{Cuadrados m\'inimos lineales}

El m\'etodo de cuadrados m\'inimos concede mayor valor relativo al 
punto que est\'a alejado del resto de los datos, pero a su vez no 
permite que domine enteramente la aproximaci\'on. 
Al carecer de la funci\'on m\'odulo (que dificulta la derivabilidad)
es una elecci\'on adecuada para trabajar. 
Adem\'as existen resultados de probabilidad y estad\'istica que 
respaldan la elecci\'on de los cuadrados m\'inimos. 

\par
En cuadrados m\'inimos lineales, se aproximan los puntos dados 
mediante una funci\'on cuyos coeficientes a determinar aparecen en 
forma lineal. 
Esto se debe a que al realizar los c\'alculos con este tipo de 
funciones, se obtiene un sistema lineal.

\par
El problema general de minimizar la suma de las diferencias al 
cuadrado en funci\'on de los coeficientes de la funci\'on $p(t)$ 
puede resolverse derivando respecto de cada uno de los coeficientes 
e igualando a cero las derivadas.

%~ , con lo que se llega a las llamadas ecuaciones normales, 
%~ sistema de ecuaciones cuya soluci\'on son los coeficientes 
%~ buscados que minimizan el error.

\newpage
\textbf{Interpretaci\'on Matricial}
\par
Se tiene un conjunto de $m$ mediciones de la forma $(x_{i}, y_{i})$ 
que se desea aproximar por una funci\'on modelo $p(t)$ que puede ser 
expresada como $\displaystyle p(t) = \sum_{j=0}^n a_j\Phi_j(t)$ donde 
$\{\Phi_0,\ldots,\Phi_n\}$ es un conjunto linealmente independiente 
de funciones conocidas.
En el caso particular que $p(t)$ sea un polinomio $\Phi_i(t) = t^i$.

\paragraph{}
Lo que se desea entonces es minimizar 
$\displaystyle f(x)=\|Ax-b\|_{2}^{2}$, donde $A$ es la matriz que 
contiene las constantes que multiplican a los coeficientes a 
determinar, $x$ tiene los coeficientes y $b$ consiste de los $y_i$. 

\paragraph{}
El sistema $Ax = b$ queda definido como
	$\begin{pmatrix}
	\Phi_0(t_1) & \Phi_1(t_1) & \cdots & \Phi_n(t_1) \\
	\Phi_0(t_2) & \Phi_1(t_2) & \cdots & \Phi_n(t_2) \\
	\vdots & \vdots & \ddots & \vdots \\
	\Phi_0(t_m) & \Phi_1(t_m) & \cdots & \Phi_n(t_m) \\
	\end{pmatrix}
	\begin{pmatrix} a_0 \\ a_1 \\ \vdots \\ a_n \end{pmatrix}
	=
	\begin{pmatrix} y_0 \\ y_1 \\ \vdots \\ y_m \end{pmatrix}$

\paragraph{}
En la mayor\'ia de los casos el sistema $Ax = b$ tal como est\'a no 
tendr\'a soluci\'on por estar \emph{sobredeterminado}. 
Por lo tanto se intenta encontrar la soluci\'on m\'as cercana, 
minimizando $f$. \\
Puede probarse que $\nabla f(x)=0$ si y s\'olo si $A^{t}Ax=A^{t}b$. \\

\textbf{Interpretaci\'on Geom\'etrica}
\par
Se busca $\displaystyle \min_x\|Ax-b\|_{2}^{2} = 
\min_{y\; \in\; Im(A)}\|y-b\|_{2}^{2}$.
Si $b\in Im(A)$ existe un $x$ tal que $Ax = b$ y por lo tanto $y - b$ 
es cero, es decir, la funci\'on modelo elegida para explicar los datos 
coincide con cada uno de ellos. 
En caso contrario, interesa buscar el $y\in Im(A)$ m\'as cercano a $b$, 
que ser\'a la proyecci\'on ortogonal de $b$ sobre la imagen de $A$, 
como puede verse en la figura.
En conclusi\'on, el $x$ buscado es el que cumple $Ax = b_{1}$ donde 
$b_{1}$ es la proyecci\'on ortogonal de $b$ en $Im(A)$.

%~ \begin{figure}[H] \centering
	%~ \subfloat{\includegraphics[scale=.20]{images/cmrecta.png}}		
	%~ \caption{Aproximaci\'on de puntos mediante una recta por 
	%~ cuadrados m\'inimos lineales.}	
%~ \end{figure}

\begin{figure}[H] \centering
	\subfloat{\includegraphics[scale=.80]{images/cmgeometrica.pdf}}		
	\caption{Interpretaci\'on geom\'etrica de Cuadrados M\'inimos.}	
\end{figure}

Para que $\|b - y\|_2$ sea m\'inima, con $y \in S$ subespacio, es 
necesario que $b - y$ pertenezca al complemento ortogonal de $S$. 
Es decir, que $y$ sea la proyecci\'on ortogonal de $b$ sobre $S$, 
en este caso la imagen de $A$. 
Por lo tanto, $b - Ax$ debe pertenecer a $Im(A)^{\bot} = Null(A^t)$. 
Para que eso suceda, debe pasar que $A^t(Ax-b)=0$. 
%~ que es justamente la soluci\'on al problema de cuadrados m\'inimos.

\subsection{Usando las ecuaciones normales}
Del enfoque anterior se desprende que el problema de cuadrados 
m\'inimos siempre tiene soluci\'on, que es \'unica si y s\'olo si 
$Null(A)=\{0\}$. 
Cuando $Ax = b$ tiene \'unica soluci\'on entonces lo mismo vale para 
$A^{t}Ax = A^{t}b$ y la soluci\'on es la misma. 
Esto permite tratar el problema de cuadrados m\'inimos con este enfoque. 

\paragraph{}
Una ventaja es que la matriz $A^{t}A$ es cuadrada, sim\'etrica y semi 
definida positiva (cuando la solucion es \'unica es definida positiva).

Sin embargo, cuando la matriz $A$ est\'a mal condicionada el uso de las 
ecuaciones normales puede empeorar el sistema, ya que $A^{t}A$ 
estar\'a peor condicionada que $A$.
Por esta raz\'on existen m\'etodos alternativos num\'ericamente m\'as 
estables.

\newpage
\subsection{Usando la factorizaci\'on QR}
Si $Q$ es una matriz ortogonal, resulta equivalente minimizar la norma 
de $Ax-b$ y la de $Q^t(Ax-b)$. 
Para aprovechar esto, se busca la descomposicion QR de la matriz $A$ 
para resolver m\'as f\'acilmente el sistema. 
Seg\'un el rango de la matriz $A$ se tienen dos casos para caracterizar 
mejor el conjunto de soluciones.
En ambos casos, se plantea el nuevo sistema $Q^t A x = Q^t b$ que es 
equivalente a $Rx = c$, con $c = (\hat{c}, d)$ donde $\hat{c}$ son 
los primeros $m$ elementos de $c$ y $d$ los restantes. 
Se busca minimizar el cuadrado del residuo $s$, que resulta ser 
$s = c - Rx$.

\subsubsection{Rango completo}
%~ En este caso $R$ es no singular. 
Los primeros $m$ elementos de $s$ son iguales a $\hat{c} - \hat{R}x$ 
y los restantes a $d$. 
 De esta forma, el cuadrado del residuo es igual a 
 $\displaystyle \|s\|^2_2 = \|\hat{c} - \hat{R}x\|^2_2 + \|d\|^2_2$.

\paragraph{}
Puesto que el t\'ermino $d$ no depende de $x$, se busca minimizar el 
primer t\'ermino. Como $\hat{R}$ era no singular, entonces la 
soluci\'on del sistema $\hat{R}x = \hat{c}$ es \'unica y es la 
soluci\'on de cuadrados m\'inimos.
%~ Cabe destacar que el t\'ermino $\|d\|^2_2$ es la norma del residuo asociado con la soluci\'on obtenida.

\subsubsection{Rango incompleto}

En este caso es necesario realizar pivoteo de columnas en la 
descomposici\'on QR. 
En cada iteraci\'on del algoritmo, se toma la columna de mayor norma, 
para dejar las columnas iguales a cero para el final. 
El algoritmo encuentra ceros luego de $r$ iteraciones, siendo $r$ el 
rango de la matriz. Resulta entonces $AP = QR$, siendo $P$ una matriz 
de permutaci\'on, $R$ una matriz con ceros debajo de la fila $r$, cuyo 
menor principal $R_r$, o $R_{11}$ es una matriz triangular superior 
no singular.

\paragraph{}
Los primeros $r$ elementos del residuo $s$ son iguales a 
$\hat{c} - R_{11} x_1 - R_{12} x_2$, siendo $x = (x_1, x_2)^t$, con 
$x_1\in\mathbb{R}^r$; y los restantes iguales a $d$. 
Nuevamente se busca minimizar el primer t\'ermino.

Se busca un $x_2$ cualquiera, a partir del cual hay un \'unico $x_1$ 
tal que $R_{11} x_1 = \hat{c} - R_{12} x_2$. 
Puesto que $R_{11}$ es no singular, y triangular superior, el sistema 
puede resolverse mediante backward substitution. \\
Se obtienen as\'i infinitas soluciones.

\begin{figure}[H] \centering
	\subfloat
	{
		$ R = \begin{pmatrix} \hat{R} \\ 0 \end{pmatrix} =
		\begin{pmatrix}
		\hat{r}_{11} & \hat{r}_{12} & \hat{r}_{13} \\
		0 & \hat{r}_{22} & \hat{r}_{23} \\
		0 & 0 & \hat{r}_{33} \\
		0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}
		c = \begin{pmatrix} \hat{c} \\ d \end{pmatrix} =
		\begin{pmatrix}
		\hat{c}_{1} \\ \hat{c}_{2} \\ \hat{c}_{3} \\
		d_{1} \\ d_{2} \\ d_{3} \end{pmatrix}$
	}
	\hspace*{30pt}
	\subfloat{$ R = \begin{pmatrix} R_{11} & R_{12} \\ 0 & 0 \end{pmatrix}$}
	\caption{Rango completo e incompleto respectivamente.}
\end{figure}

Usar la descomposici\'on QR es poco estable, y hallar el rango de $A$ 
no siempre es sencillo por errores de redondeo. 
Por este motivo existe otra descomposici\'on, 
Singular Value Descomposition (SVD).

\subsection{Usando la descomposici\'on SVD}

\begin{theorem}
	Sea $A$ en $\real^{m\times n}$ con $rango(A)=r$, entonces existen 
	$\{v_1, \ldots, v_n\}$ y $\{u_1,\ldots,u_m\}$ tales que forman 
	una base ortonormal de $\real^n$ y de $\real^m$ respectivamente, 
	de forma que $U^t A V = D$.\\
	$D = \begin{pmatrix} D_r & 0 \\ 0 & 0 \end{pmatrix}$
	donde $D_r = \textrm{diag}(s_1,\ldots,s_r)$ y $s_i$ son los valores 
	singulares de $A$ ordenados de mayor a menor. Los $u_i$ son los 
	autovectores normalizados de $AA^t$, y los $v_i$ son los 
	autovectores normalizados de $A^tA$.
	
	\paragraph{}
	Esto implica que $A = UDV^t$, o que $AV = UD$, siendo $U$ y $V$ 
	matrices ortogonales de $m\times m$ y $n\times n$ respectivamente 
	y $D$ una matriz diagonal de $m\times n$, misma dimensi\'on que $A$. 
	Esta descomposici\'on existe para cualquier matriz $A$.
\end{theorem}

\newpage
Con el mismo criterio que en $QR$, puesto que $U$ es ortogonal, 
minimizar la norma de $b - Ax$ es equivalentente a minimizar la norma 
de $U^t(b-Ax) = U^tb - DV^tx$. 
Tomando $c = U^tb$ e $y = V^tx$, minimizar la norma de $b - Ax$ 
equivale a minimizar la de $c - Dy$. \\
Calculando, $\displaystyle \|b-Ax\|^2_2 = \|c-Dy\|^2_2 = 
\sum_{i=1}^r{|c_i - s_iy_i|^2} + \sum_{i=r+1}^m{|c_i|^2}. $

Por lo tanto, la soluci\'on del sistema es $\displaystyle y_i = 
\frac{c_i}{s_i} \; \forall \; i = 1,\ldots, r$.
Si $r < m$, la matriz original $A$ no tenia rango completo, entonces 
los valores $y_i$ para $i = r+1 \ldots m$ no est\'an involucrados en 
la expresi\'on a minimizar, y pueden tomar cualquier valor. 
Si se busca el $x$ de norma m\'inima resultante de $x = Vy$, entonces 
se toman los $y_i$ restantes iguales a cero.

\paragraph{}
El m\'etodo SVD es m\'as computacionalmente costoso que usando QR, 
independientemente de c\'omo se organicen los c\'alculos. 
Se pueden hacer optimizaciones, como no calcular nunca la matriz $U$ 
sino aplicar solamente los reflectores sobre $b$, o calcular solamente 
las primeras $r$ columnas de $V$ que ser\'an las necesarias para 
resolver $x = Vy$, con los $y_i$ $i>r$ iguales a cero. 
De todas formas, la estabilidad de SVD en casos de rango incompleto 
lo hace preferible por sobre QR.

\begin{theorem}
	La estabilidad de la solucion de cuadrados minimos lineales esta 
	dada por la siguiente cota, siendo $\bar{x}$ la soluci\'on al 
	problema de cuadrados m\'inimos y $\bar{b}$ el proyector sobre 
	la imagen de $A$, 
	
	\begin{displaymath}
		\varepsilon_r (\bar{x}) \leq \chi(A) \varepsilon_r (\bar{b})	
	\end{displaymath}
	
	Donde $\chi(A) = \|A\|_2 \|(A^tA)^{-1}\|_2$ es una generalizaci\'on 
	del n\'umero de condici\'on, cuando $A$ no tiene necesariamente 
	inversa.
\end{theorem}
